Analysis_OleateMicrobiota

Project Description

In this current project, we have 3 groups of mice. A control group that was fed 10% fat diet, an oleate group which was fed with an oleate diet and a HFD group which was fed with 60% fat diet. The project included 19 mice. We collected feces, contents from the ileum, cecum and colon and we performed 16S rRNA sequencing on them. The DNA extraction and sequencing was done in the university of Michigan Microbiome core. The details of sequencing procedures are present in the Microbiome_core_sequencing_details.xlsx excel file.

Preprocessing and classifying OTUs in mothur

First, we created a file that assign the fastq files to the samples. The file is called oleate.files in the data directory.
We used mothur to preprocess the data and to clasify OTUs. An executable version of mothur for windows (version 1.44.3) was installed in the folder containing the fastq files. The original fastq files were processed using mothur pipeline using the following code:

# setting the input and output directories  
set.dir(input=Data, output= ButyrateResults)
# making contigs  
make.contigs(file=oleate.files, processors=4)
screen.seqs(fasta=current, group=current, maxambig=0, maxlength=275)
unique.seqs(fasta=current)
count.seqs(name=current, group=current)
summary.seqs(count=current)
# make aan align reference file from silva original reference file. The silva reference files (version 138.1) were downloaded from <https://mothur.org/wiki/silva_reference_files/>
pcr.seqs(fasta=silva.nr_v138_1.align, start=11894, end=25319, keepdots=F, processors=8)
# align sequences  
align.seqs(fasta= oleate.trim.contigs.good.unique.fasta, reference= silva.nr_v138_1.pcr.align)
summary.seqs(fasta= current, count= current)
screen.seqs(fasta=current, count=current, summary=current, start=1968, end=11550, maxhomop=8)
summary.seqs(fasta= current, count= current)
filter.seqs(fasta=current, vertical=T, trump=.)
unique.seqs(fasta=current, count=current)
pre.cluster(fasta=current, count=current, diffs=2)
chimera.vsearch(fasta=current, count=current, dereplicate=t)
remove.seqs(fasta=current, accnos=current)
summary.seqs(fasta=current, count=current)
# classify sequences
classify.seqs(fasta=current, count=current, reference= silva.nr_v138_1.pcr.align, taxonomy= silva.nr_v138_1.tax, cutoff=80)
remove.lineage(fasta=current, count=current, taxonomy=current, taxon=Chloroplast-Mitochondria-unknown-Archaea-Eukaryota)
summary.tax(taxonomy=current, count=current)
rename.file(taxonomy=oleate.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.0.03.cons.taxonomy, shared= oleate.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.opti_mcc.shared)
# investigate the read counts in each sample and selecting the least count for subsampling
count.groups(shared=oleate.opti_mcc.shared)
sub.sample(shared= current, size= 2766)
rarefaction.single(shared=current, calc=sobs, freq=100)
summary.single(shared=current, calc=nseqs-coverage-sobs-invsimpson, subsample=T)
# calculating distances and principal coordinate analysis
dist.shared(shared=current, calc=thetayc-jclass, subsample=t)
pcoa(phylip=current)
nmds(phylip=current)
nmds(phylip=current, mindim=3, maxdim=3)
# we created a set of txt files in excel to reflect the design and metadata of the samples. these files are present in the data directory under the names `condition.design`, `site.design` and `oleate.metadata`. Next, we proceeded with performing AMOVA to detect significance
amova(phylip=current, design=condition.design)
amova(phylip=current, design=site.design)
homova(phylip=current, design=condition.design)
homova(phylip=current, design=site.design)

Further analysis and visualization in R

We created an R project and copied the following files to the project directory:
1- oleate.opti_mcc.shared (the OTU table of all samples)
2- oleate.opti_mcc.0.03.subsample.shared (the OTU table after subsampling)
3-oleate.metadata (a txt file that includes the sample metadata)
4- oleate.taxonomy (the taxonomy file that was generated by mothur. we modified this file in excel. As it is from mothur, there are not column headers for order, family genus etc. in the cons.taxonomy file and this must be fixed first. This is how I did it: read cons.taxonomy file into excel and choose semicolon as a separator so each tax level is therefore in its own column. then delete header ‘taxonomy’ and put in appropriate header for each column (order or family or genus etc. Then select all, press ctrl+shift+H on mac to find and replace. Replate all (*) with nothing)
## Install and load required libraries

library("phyloseq")
library("ggplot2")
library("plyr")
library("vegan")
library("grid")
library("directlabels")
library("knitr")
library("clustsig")
library("ellipse")
library("ape")
library("DESeq2")
library("microbial")
library("VennDiagram")
library("microbiome")
library("microbiomeutilities")

## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'

library("viridis")
library("RColorBrewer")
library("ggpubr")
library("patchwork")

Make a phyloseq object from mothur files

## import data 
sharedFile = read.table("oleate.opti_mcc.shared")
sharedFile = t(sharedFile) ## transform data
rownames(sharedFile) = sharedFile[,1]  ## define rowname
colnames(sharedFile) = sharedFile[2,]  ## define column names
sharedFile <-as.matrix(sharedFile)
sharedFile = sharedFile[,2:57]       ## only include columns that reflects OTU numbers (remember order [row,colum]) 2 -> 57 for me
## as I had 57 samples
sharedFile = sharedFile[4:819,]     ## and which rows you want. I have 819 OTUs so I changed this number to 819
# original => sharedFile = sharedFile[4:819,]
class(sharedFile) <- "numeric"
head(sharedFile)  ## look at head first 7 lines to see it's made correclty

##        ce4406 ce4407 ce4408 ce4410 ce4411 ce4412 ce4413 ce4414 ce4415 ce4416
## Otu001    648    490    686    411    726    838    761    175    144   1520
## Otu002      1    143   1693   1037    794    548    356    418    149   1075
## Otu003    867    743   1670    733   1171    559    953    165    202    564
## Otu004    310      8   1379    808    903    419   1487    953     46    192
## Otu005    483    704   1357   2703    825    735    381    416    188   2496
## Otu006      0      0      0      1    170      2      4     49     40    365
##        ce4417 ce4418 ce4419 ce4420 ce4421 ce4422 ce4423 ce4424 co4406 co4407
## Otu001   1118   1091   1691    703    936   1011     54   1452   1170    257
## Otu002   2412   6002   1290    586   1216   1117   2346     48      5    843
## Otu003   1201   1204    871    980    890    922   1522    388    406   1539
## Otu004     29    291   1469     21      6      4      3     20   1066     14
## Otu005   2134    969   1805   1737   1385   1602   1220    391    292    188
## Otu006    123   1317     97    245    237    272    146   2061      3      2
##        co4408 co4410 co4411 co4412 co4413 co4414 co4415 co4416 co4417 co4418
## Otu001   1823    856    397    705    356    392     45    201   4331    356
## Otu002    767    177    220    404     24    442    261   1187   1901   1269
## Otu003   1116    517    909    996    394    499    490    498    664   1311
## Otu004   1089    564    518    468    660   1151     27    133    481     13
## Otu005    264    400    795    778    310    875    168     31    804    483
## Otu006      0      0     59      1      1    110     19    195   1149     66
##        co4419 co4420 co4421 co4422 co4423 co4424 co4428 f4406 f4407 f4408 f4410
## Otu001    182    188   1372    944   1621    233   1166  7021  9127  2210  3181
## Otu002   2380   1814    815    390    778   2231     10    31  1790    14   173
## Otu003    567    499    527    791   1318   1249    225  1108  1447   679   726
## Otu004     31    484     35      0     18      4      6  1221  2604   386  2143
## Otu005    281    212    501    519    584    237    148  1300   380  1279   628
## Otu006    240     39    280    237    667    376   1572     3     0     0     1
##        f4411 f4412 f4413 f4414 f4415 f4416 f4417 f4418 f4419 f4420 f4421 f4422
## Otu001  6323  8943  2011  4584  4530  7401  3309  6001  7938  4862  1619  1334
## Otu002  1367    30    99  2063  1193   898   472   265  3751  3667  1244  1022
## Otu003  1242  1133  1532  1163   687   902   536   443  1878   554  1781  1668
## Otu004  2441   716  1082  4784  4027  3709  3814  2177  1569  2593    80    21
## Otu005  1154  1847  1024   816   591  1263  1511   521  1505  1013   272   558
## Otu006   786     1     0  1648   862  1234  1175   931  2214  1074  2252  1979
##        f4423 f4424 f4428
## Otu001  3749 11842   791
## Otu002  1421  1194   996
## Otu003   610  2459   971
## Otu004   147    23     5
## Otu005  1367    64   324
## Otu006   798   506  1765

## Import subsampled otu matrix (26880 seqs)
sharedsubFile = read.table('oleate.opti_mcc.0.03.subsample.shared')
sharedsubFile = t(sharedsubFile)
rownames(sharedsubFile) = sharedsubFile[,1]
colnames(sharedsubFile) = sharedsubFile[2,]
sharedsubFile = sharedsubFile[,2:57]
sharedsubFile = sharedsubFile[4:494,]
class(sharedsubFile) <- "numeric"
head(sharedsubFile)

##        ce4406 ce4407 ce4408 ce4410 ce4411 ce4412 ce4413 ce4414 ce4415 ce4416
## Otu001    335    271    174     72    204    358    233    131    144    312
## Otu002      1     68    361    199    235    236    103    311    149    237
## Otu003    468    380    379    139    359    254    289    125    202    108
## Otu004    145      2    312    128    261    186    453    680     46     37
## Otu005    245    348    299    492    254    319    124    308    188    509
## Otu006      0      0      0      0     38      1      1     31     40     84
##        ce4417 ce4418 ce4419 ce4420 ce4421 ce4422 ce4423 ce4424 co4406 co4407
## Otu001    187    161    311    109    212    197     10    137    456    145
## Otu002    414    825    219     86    255    194    590      3      2    528
## Otu003    213    166    156    157    197    176    380     30    160    916
## Otu004      6     38    284      3      0      1      0      1    407      9
## Otu005    412    134    311    262    300    295    299     34    128    109
## Otu006     20    173     20     46     47     44     30    184      2      2
##        co4408 co4410 co4411 co4412 co4413 co4414 co4415 co4416 co4417 co4418
## Otu001    473    406    172    289    160    168     41    127    611    138
## Otu002    224     82    101    156     11    181    240    770    281    465
## Otu003    304    272    445    419    180    214    459    308     85    543
## Otu004    317    258    253    167    321    469     27     86     69      7
## Otu005     83    202    400    294    147    373    155     19    105    172
## Otu006      0      0     24      0      0     44     18    121    169     18
##        co4419 co4420 co4421 co4422 co4423 co4424 co4428 f4406 f4407 f4408 f4410
## Otu001     84     80    358    335    340     72    237   933  1161   548   615
## Otu002    873    792    208    118    166    744      1     4   230     2    29
## Otu003    212    202    130    290    274    413     48   141   199   168   115
## Otu004     18    223      9      0      4      1      2   151   360   105   410
## Otu005    111     73    132    184    122     82     24   154    52   311   113
## Otu006     95     14     68     96    126    127    293     0     0     0     0
##        f4411 f4412 f4413 f4414 f4415 f4416 f4417 f4418 f4419 f4420 f4421 f4422
## Otu001   929  1220   419   585   727   779   584  1196   788   533   250   228
## Otu002   196     5    22   259   167   105    74    49   390   416   224   144
## Otu003   171   167   296   144   124    95    99    87   195    75   271   261
## Otu004   384    99   204   639   655   415   673   421   145   281    15     5
## Otu005   179   283   209    99   101   125   256    93   145   131    40    87
## Otu006   106     0     0   177   135   128   212   197   229   116   332   338
##        f4423 f4424 f4428
## Otu001   401  1571   135
## Otu002   158   159   161
## Otu003    65   348   150
## Otu004    12     3     0
## Otu005   137    11    59
## Otu006    95    66   275

dim(sharedsubFile)

## [1] 491  56

sharedsubFile <-as.matrix(sharedsubFile)
## Import taxonomy file 
## Copy and paste into notepad and save. Then carry on: 
taxFile = read.table('oleate.taxonomy', header=T, sep='\t')
rownames(taxFile) = taxFile[,1]
taxFile = taxFile[,3:8]
taxFile = as.matrix(taxFile)
head(taxFile)

##        Kingdom    Phylum              Class              Order               
## Otu001 "Bacteria" "Firmicutes"        "Bacilli"          "Lactobacillales"   
## Otu002 "Bacteria" "Verrucomicrobiota" "Verrucomicrobiae" "Verrucomicrobiales"
## Otu003 "Bacteria" "Bacteroidota"      "Bacteroidia"      "Bacteroidales"     
## Otu004 "Bacteria" "Firmicutes"        "Bacilli"          "Lactobacillales"   
## Otu005 "Bacteria" "Firmicutes"        "Clostridia"       "Lachnospirales"    
## Otu006 "Bacteria" "Firmicutes"        "Bacilli"          "Erysipelotrichales"
##        Family                Genus            
## Otu001 "Streptococcaceae"    "Lactococcus"    
## Otu002 "Akkermansiaceae"     "Akkermansia"    
## Otu003 "Bacteroidaceae"      "Bacteroides"    
## Otu004 "Lactobacillaceae"    "Lactobacillus"  
## Otu005 "Lachnospiraceae"     "Lachnospiraceae"
## Otu006 "Erysipelotrichaceae" "Dubosiella"

## import metadata file
metaFile = read.table('oleate.metadata', header=T, sep='\t')
rownames(metaFile) = metaFile[,1]
metaFile = metaFile[,1:4]
head(metaFile)

##         group  site mouse condition
## ce4406 ce4406 cecum  4406    10%Fat
## ce4407 ce4407 cecum  4407    10%Fat
## ce4408 ce4408 cecum  4408    10%Fat
## ce4410 ce4410 cecum  4410    10%Fat
## ce4411 ce4411 cecum  4411    10%Fat
## ce4412 ce4412 cecum  4412    10%Fat

## Create phyloseq object
OTU = otu_table(sharedFile, taxa_are_rows = TRUE)
OTUsub = otu_table(sharedsubFile, taxa_are_rows = TRUE)
TAX = tax_table(taxFile)
META = sample_data(metaFile)
physeq = phyloseq(OTU, TAX, META)
physeqSub = phyloseq(OTUsub, TAX, META)
physeqSub

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 491 taxa and 56 samples ]
## sample_data() Sample Data:       [ 56 samples by 4 sample variables ]
## tax_table()   Taxonomy Table:    [ 491 taxa by 6 taxonomic ranks ]

## Get rid of any OTUs not present in any samples and get relative abundance
microSub <- prune_taxa(taxa_sums(physeqSub) > 0, physeqSub)
microSubRel = transform_sample_counts(microSub, function(x) x / sum(x) )
microSubRelFilt = filter_taxa(microSubRel, function(x) mean(x) > 1e-5, TRUE)
# for subsampled shared file
# create a tree and a new phyloseq object
random_tree = rtree(ntaxa(microSub), rooted=TRUE, tip.label=taxa_names(microSub))
plot(random_tree)

microSubt = merge_phyloseq(microSub, random_tree)

Plotting core graphs

# Plotting rarefaction
rarecurve(t(otu_table(microSubt)), step=50, cex=0.5)

# Plot a prettier rarefaction curve
# set seed
set.seed(1)
subsamples <- seq(0, 5000, by=100)[-1]
#subsamples = c(10, 5000, 10000, 20000, 30000)
p <- plot_alpha_rcurve(microSub, index="observed", 
                       subsamples = subsamples,
                       lower.conf = 0.025, 
                       upper.conf = 0.975,
                       group="condition",
                       label.color = "brown3",
                       label.size = 3,
                       label.min = TRUE) 
# change color of line 
mycols <- c("brown3", "steelblue","grey50")

p <- p + scale_color_manual(values = mycols) + 
  scale_fill_manual(values = mycols)
print(p)

# Plotting taxonomy
plot_bar(microSubt, fill="Phylum") + 
  facet_wrap(~condition, scales = "free_x", nrow = 1) +
  theme(text = element_text(size=16))+geom_bar(stat="identity") +
  scale_fill_brewer(type = "seq", palette = "Spectral") +
  theme(
    legend.background = element_rect(
      fill = "lemonchiffon", 
      colour = "grey50", 
      size = 1
    )
  )

plot_bar(microSubt, fill="Phylum") + 
  facet_wrap(~site, scales = "free_x", nrow = 1) +
  theme(text = element_text(size=16))+geom_bar(stat="identity") +
  scale_fill_brewer(type = "seq", palette = "Spectral") +
  theme(
    legend.background = element_rect(
      fill = "lemonchiffon", 
      colour = "grey50", 
      size = 1
    )
  ) 

plot_bar(microSub, fill="Phylum") +
  facet_grid(condition~site, scales = "free") +  ## separate facets
  theme(text = element_text(size=16)) +  ## increase font of all elements
  geom_bar(stat="identity")+
  scale_fill_brewer(type = "seq", palette = "Spectral") +
  theme(
    legend.background = element_rect(
      fill = "lemonchiffon", 
      colour = "grey50", 
      size = 1
    )
  )

# Plotting PCoA
my.PCoA2 <- ordinate(microSub, "PCoA", "bray")
plot_ordination(microSub, my.PCoA2, type = "group", color = "condition", shape = "site")+theme(text = element_text(size=20))+ geom_point(size=5) + scale_color_manual(values=c("brown3", "steelblue","grey50"))

plot_ordination(microSub, my.PCoA2, type = "group", color = "condition")+ facet_wrap(~site, scales = "free_x", nrow = 1) + geom_point(size=5)+theme(text = element_text(size=20)) + scale_color_manual(values=c("brown3", "steelblue","grey50"))

plot_ordination(microSub, my.PCoA2, type = "group", color = "condition", shape = "site")+theme(text = element_text(size=20))+ geom_point(size=5)+ 
  geom_line() + scale_color_manual(values=c("brown3", "steelblue","grey50"))

plot_ordination(microSub, my.PCoA2, type = "group", color = "condition")+theme(text = element_text(size=20))+ geom_point(size=5)+ stat_ellipse(level=0.9)  + scale_color_manual(values=c("brown3", "steelblue","grey50"))

# PCoA plot using the unweighted UniFrac as distance
wunifrac_dist = phyloseq::distance(microSubt, method="unifrac", weighted=F)
ordination = ordinate(microSubt, method="PCoA", distance=wunifrac_dist)
plot_ordination(microSubt, ordination, color="condition") + theme(aspect.ratio=1)+
  ggtitle("PCoA: unweigthed Unifrac")+geom_point(size=5)  + scale_color_manual(values=c("brown3", "steelblue","grey50"))

#Plot PCoA using the weighted UniFrac as distance
wunifrac_dist = phyloseq::distance(microSubt, method="unifrac", weighted=T)
ordUF = ordinate(microSubt, method="PCoA", distance=wunifrac_dist)
plot_ordination(microSubt, ordUF, color = "condition") + 
  ggtitle("PCoA: weigthed Unifrac")+geom_point(size=5) + scale_color_manual(values=c("brown3", "steelblue","grey50"))

## Normalization and plotting figures by group

# Normalization
phy <- normalize(microSubt, method = "relative")

## Normalization using relative method

# Plot taxonomy by group
plotbar(phy,level="Phylum", group="condition") + theme(axis.text.x = element_text(size = 16))

plotbar(phy,level="Phylum", group="site") + theme(axis.text.x = element_text(size = 16))

# Plot alpha diversity
plotalpha(microSubt, group = "condition", color = c("brown3","grey50", "steelblue"))

plotalpha(microSubt, group = "site")

# Test for significance between groups
beta_condition <-betatest(phy,group="condition")

## Do PERMANOVA for: condition

beta_site <-betatest(phy,group="site")

## Do PERMANOVA for: site

# Biomarkers discovery and plotting
bio <- biomarker(microSubt,group="condition")

## Normalization using relative method

plotmarker(bio,level="Genus")

# LefSe testing and plotting
lda <- ldamarker(microSubt,group="condition")

## Normalization using relative method

write.csv(lda, "lda.csv", row.names = FALSE)
# open the csv file in excel, separate the taxonomy to kingdom, phylum,... and save it as lda2.csv so that the names can appear better on the graph
lda2 <- read.csv(file = 'lda2.csv')
plotLDA(lda2,group=c("10%Fat","60%Fat"),lda=5,pvalue=0.05)+theme(axis.text.y = element_text(size = 16))+theme(text = element_text(size=20))

plotLDA(lda2,group=c("10%Fat","Oleate"),lda=5,pvalue=0.05)+theme(axis.text.y = element_text(size = 16))+theme(text = element_text(size=20))

plotLDA(lda2,group=c("Oleate","60%Fat"),lda=5,pvalue=0.05) +theme(axis.text.y = element_text(size = 16))+theme(text = element_text(size=20))

## Modified difftest function from microbial package to test for significance

# create difftest2 function
difftest2<-function(physeq,group,pvalue=0.05,padj=NULL,log2FC=0,gm_mean=TRUE,fitType="local",quiet=FALSE){
  if(!taxa_are_rows(physeq)){
    physeq<-t(physeq)
  }
  otu <- as(otu_table(physeq),"matrix")
  tax <- as.data.frame(as.matrix(tax_table(physeq)))
  colData<-as(sample_data(physeq),"data.frame")
  colData$condition<-colData[,group]
  contrasts<-levels(factor(unique(colData$condition)))[1:2]
  if(isTRUE(gm_mean)){
    countData<-round(otu, digits = 0)
  }else{
    countData<-round(otu, digits = 0)+1
  }
  dds <- DESeqDataSetFromMatrix(countData, colData, as.formula(~condition))
  if(isTRUE(gm_mean)){
    geoMeans = apply(counts(dds), 1, gm_mean)
    dds = estimateSizeFactors(dds, geoMeans = geoMeans)
  }
  dds <- DESeq(dds, fitType=fitType)
  res <- results(dds,contrast=c("condition",contrasts),cooksCutoff = FALSE)
  res_tax = cbind(as.data.frame(res), as.matrix(countData[rownames(res), ]))
  if(!is.null(padj)){
    pval<-padj
    sig <- rownames(subset(res,padj<pval&abs(log2FoldChange)>log2FC))
  }else{
    pval<-pvalue
    sig <- rownames(subset(res,pvalue<pval&abs(log2FoldChange)>log2FC))
  }
  res_tax$Significant<- "No"
  res_tax$Significant <- ifelse(rownames(res_tax) %in% sig, "Yes", "No")
  res_tax <- cbind(res_tax, tax[rownames(res),])
  return(as.data.frame(res_tax))
}
# Test for significant bacteria using microbial package
mic_res <- difftest2(microSubt,group="condition", gm_mean=FALSE)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## estimating size factors

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## final dispersion estimates

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## fitting model and testing

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## fitting model and testing

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

plotdiff(mic_res,level="Genus",padj=0.001,log2FC = 3,fontsize.y = 14)

## Test for significant OTUs using DESeq2

sample_data(microSubt)$condition <- as.factor(sample_data(microSubt)$condition)
ds = phyloseq_to_deseq2(microSubt, ~ condition)

## converting counts to integer mode

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

ds = DESeq(ds)

## estimating size factors
##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## final dispersion estimates

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## fitting model and testing

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## fitting model and testing

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

alpha = 0.01
# compare Oleate and 60%Fat
res = results(ds, contrast=c("condition", "Oleate", "60%Fat"), alpha=alpha)
res = res[order(res$padj, na.last=NA), ]
res_sig = res[(res$padj < alpha), ]
res_sig

## log2 fold change (MLE): condition Oleate vs 60%Fat 
## Wald test p-value: condition Oleate vs 60%Fat 
## DataFrame with 25 rows and 6 columns
##         baseMean log2FoldChange     lfcSE      stat      pvalue        padj
##        <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
## Otu004 200.56502        6.13780  0.634629   9.67149 3.98541e-22 7.81140e-20
## Otu031  17.90550        3.11457  0.399384   7.79842 6.26851e-15 6.14314e-13
## Otu028  18.25563        4.51975  0.610052   7.40880 1.27449e-13 8.32670e-12
## Otu029  21.61460       -5.03417  0.761618  -6.60984 3.84734e-11 1.80406e-09
## Otu054   7.35816        2.95435  0.448766   6.58327 4.60220e-11 1.80406e-09
## ...          ...            ...       ...       ...         ...         ...
## Otu050   8.71757       -2.90784  0.699150  -4.15910 3.19497e-05 0.000298197
## Otu027  21.35176        2.17939  0.534261   4.07927 4.51775e-05 0.000402491
## Otu040  11.60680        1.84137  0.487587   3.77648 1.59059e-04 0.001355458
## Otu078   2.00088       -3.40305  0.941310  -3.61523 3.00082e-04 0.002450666
## Otu062   5.35700        1.88742  0.584235   3.23058 1.23540e-03 0.009685536

res_sig = cbind(as(res_sig, "data.frame"), as(tax_table(microSub)[rownames(res_sig), ], "matrix"))
ggplot(res_sig, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
  geom_jitter(size=5, width = 0.2) + theme(text = element_text(size=16)) +
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))

# compare 10%Fat and 60%Fat
res2 = results(ds, contrast=c("condition", "10%Fat", "60%Fat"), alpha=alpha)
res2 = res2[order(res2$padj, na.last=NA), ]
res_sig2 = res2[(res2$padj < alpha), ]
res_sig2

## log2 fold change (MLE): condition 10%Fat vs 60%Fat 
## Wald test p-value: condition 10%Fat vs 60%Fat 
## DataFrame with 53 rows and 6 columns
##         baseMean log2FoldChange     lfcSE      stat      pvalue        padj
##        <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
## Otu007   92.6627       -8.96186  0.564355 -15.87981 8.74359e-57 1.30280e-54
## Otu008   60.6831       -5.59715  0.559648 -10.00120 1.50568e-23 7.47822e-22
## Otu036   13.0374       -7.57488  0.756046 -10.01907 1.25680e-23 7.47822e-22
## Otu033   12.2718       -7.32904  0.771735  -9.49683 2.16382e-21 8.06022e-20
## Otu029   21.6146       -7.82794  0.852785  -9.17927 4.34030e-20 1.12102e-18
## ...          ...            ...       ...       ...         ...         ...
## Otu020  28.00857        1.19042  0.404478   2.94311  0.00324931  0.00949308
## Otu018  30.12132        1.71064  0.576419   2.96771  0.00300032  0.00949308
## Otu050   8.71757       -2.02837  0.687046  -2.95231  0.00315407  0.00949308
## Otu013  64.09575        1.45752  0.497800   2.92791  0.00341245  0.00977797
## Otu012  43.93690       -1.27908  0.438006  -2.92024  0.00349767  0.00983306

res_sig2 = cbind(as(res_sig2, "data.frame"), as(tax_table(microSub)[rownames(res_sig2), ], "matrix"))
ggplot(res_sig2, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
  geom_jitter(size=5, width = 0.2) + theme(text = element_text(size=16)) +
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))

# compare 10%Fat and Oleate
res3 = results(ds, contrast=c("condition", "10%Fat", "Oleate"), alpha=alpha)
res3 = res3[order(res3$padj, na.last=NA), ]
res_sig3 = res3[(res3$padj < alpha), ]
res_sig3

## log2 fold change (MLE): condition 10%Fat vs Oleate 
## Wald test p-value: condition 10%Fat vs Oleate 
## DataFrame with 34 rows and 6 columns
##         baseMean log2FoldChange     lfcSE      stat      pvalue        padj
##        <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
## Otu007   92.6627       -9.43646  0.542183 -17.40455 7.62021e-68 1.20399e-65
## Otu035   14.4848       -7.15821  0.783378  -9.13761 6.38477e-20 5.04397e-18
## Otu036   13.0374       -6.40371  0.707215  -9.05484 1.36771e-19 7.20329e-18
## Otu033   12.2718       -6.20581  0.720641  -8.61152 7.21001e-18 2.84795e-16
## Otu023   20.7554       -4.48379  0.549579  -8.15858 3.38974e-16 1.07116e-14
## ...          ...            ...       ...       ...         ...         ...
## Otu029 21.614602       -2.79377  0.810938  -3.44510 0.000570838  0.00300641
## Otu030 16.043388        1.23116  0.375613   3.27775 0.001046392  0.00516656
## Otu068  2.640345       -3.05941  0.933212  -3.27836 0.001044109  0.00516656
## Otu125  0.721507       -2.90713  0.906565  -3.20676 0.001342406  0.00642728
## Otu019 21.243604       -2.03839  0.663115  -3.07396 0.002112387  0.00981639

res_sig3 = cbind(as(res_sig3, "data.frame"), as(tax_table(microSubt)[rownames(res_sig3), ], "matrix"))
ggplot(res_sig3, aes(x=Genus, y=log2FoldChange, color=Phylum)) +
  geom_jitter(size=5, width = 0.2) + theme(text = element_text(size=16)) +
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))

## Plot heatmap

plot_taxa_heatmap(microSubt,
                  subset.top = 20,
                  VariableA = "condition",
                  heatcolors = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(100),
                  transformation = "log10"
)

## Top 20 OTUs selected

## log10, if zeros in data then log10(1+x) will be used

## First top taxa were selected and 
## then abundances tranformed to log10(1+X)

## Warning in transform(phyobj1, "log10"): OTU table contains zeroes. Using log10(1
## + x) transform.

## Venn Diagram

pseq<-microSubt
table(meta(pseq)$condition, useNA = "always")

## 
## 10%Fat 60%Fat Oleate   <NA> 
##     21     14     21      0

# convert to relative abundances
pseq.rel <- microbiome::transform(pseq, "compositional")
disease_states <- unique(as.character(meta(pseq.rel)$condition))
print(disease_states)

## [1] "10%Fat" "Oleate" "60%Fat"

#Write a for loop to go through each of the disease_states one by one and combine identified core taxa into a list.
list_core <- c() # an empty object to store information

for (n in disease_states){ # for each variable n in DiseaseState
  #print(paste0("Identifying Core Taxa for ", n))
  
  ps.sub <- subset_samples(pseq.rel, condition == n) # Choose sample from DiseaseState by n
  
  core_m <- core_members(ps.sub, # ps.sub is phyloseq selected with only samples from g 
                         detection = 0.001, # 0.001 in atleast 90% samples 
                         prevalence = 0.75)
  print(paste0("No. of core taxa in ", n, " : ", length(core_m))) # print core taxa identified in each DiseaseState.
  list_core[[n]] <- core_m # add to a list core taxa for each group.
  #print(list_core)
}

## [1] "No. of core taxa in 10%Fat : 38"
## [1] "No. of core taxa in Oleate : 39"
## [1] "No. of core taxa in 60%Fat : 38"

print(list_core)

## $`10%Fat`
##  [1] "Otu046" "Otu015" "Otu011" "Otu037" "Otu030" "Otu045" "Otu005" "Otu001"
##  [9] "Otu012" "Otu052" "Otu042" "Otu034" "Otu028" "Otu027" "Otu014" "Otu061"
## [17] "Otu002" "Otu075" "Otu013" "Otu017" "Otu016" "Otu025" "Otu003" "Otu004"
## [25] "Otu040" "Otu020" "Otu026" "Otu048" "Otu059" "Otu024" "Otu018" "Otu031"
## [33] "Otu009" "Otu021" "Otu022" "Otu053" "Otu051" "Otu054"
## 
## $Oleate
##  [1] "Otu046" "Otu015" "Otu011" "Otu037" "Otu074" "Otu045" "Otu005" "Otu001"
##  [9] "Otu012" "Otu052" "Otu042" "Otu034" "Otu028" "Otu027" "Otu023" "Otu014"
## [17] "Otu002" "Otu013" "Otu062" "Otu039" "Otu016" "Otu025" "Otu003" "Otu004"
## [25] "Otu008" "Otu044" "Otu041" "Otu040" "Otu020" "Otu026" "Otu006" "Otu059"
## [33] "Otu024" "Otu031" "Otu009" "Otu055" "Otu021" "Otu054" "Otu007"
## 
## $`60%Fat`
##  [1] "Otu015" "Otu011" "Otu010" "Otu037" "Otu030" "Otu045" "Otu005" "Otu001"
##  [9] "Otu012" "Otu033" "Otu029" "Otu034" "Otu023" "Otu014" "Otu076" "Otu002"
## [17] "Otu013" "Otu017" "Otu035" "Otu039" "Otu025" "Otu003" "Otu008" "Otu044"
## [25] "Otu041" "Otu020" "Otu026" "Otu036" "Otu006" "Otu019" "Otu048" "Otu024"
## [33] "Otu018" "Otu078" "Otu009" "Otu021" "Otu050" "Otu007"

# Specify colors and plot venn
# supplying colors in the order they appear in list_core
mycols <- c("brown3", "grey50", "steelblue")
venn.diagram(list_core,fill = mycols, filename = '#14_venn_diagramm.png',
             output=TRUE)

## [1] 1

Plotting top 4 genus

pa <- aggregate_taxa(microSubt, "Genus")
top_four <- top_taxa(pa, 4)
top_four

## [1] "Lachnospiraceae" "Lactococcus"     "Muribaculaceae"  "Akkermansia"

mycols <- c("brown3", "steelblue", "grey50")
p <- plot_listed_taxa(pa, top_four, 
                      group= "condition",
                      group.order = c("10%Fat","Oleate","60%Fat"),
                      group.colors = mycols,
                      add.violin = TRUE,
                      violin.opacity = 0.3,
                      dot.opacity = 0.25,
                      box.opacity = 0.25,
                      panel.arrange= "wrap")

## An additonal column Sam_rep with sample names is created for reference purpose

comps <- make_pairs(sample_data(pa)$condition)
p <- p + stat_compare_means(
  comparisons = comps,
  label = "p.format",
  tip.length = 0.05,
  method = "wilcox.test")
print(p + ylab("Relative abundance") + scale_y_continuous(labels = scales::percent))

## Get taxa summary by group(s)

p0 <- microSubt
p0.gen <- aggregate_taxa(microSubt,"Genus")
x.d <- dominant_taxa(p0,level = "Genus", group="condition")
head(x.d$dominant_overview)

## # A tibble: 6 x 5
## # Groups:   condition [3]
##   condition dominant_taxa       n rel.freq rel.freq.pct
##   <fct>     <chr>           <int>    <dbl> <chr>       
## 1 10%Fat    Lachnospiraceae    14     66.7 67%         
## 2 60%Fat    Lachnospiraceae     9     64.3 64%         
## 3 Oleate    Lachnospiraceae     9     42.9 43%         
## 4 10%Fat    Lactococcus         6     28.6 29%         
## 5 Oleate    Lactococcus         6     28.6 29%         
## 6 Oleate    Akkermansia         4     19   19%

tx.sum1 <- taxa_summary(p0, "Phylum")

## Data provided is not compositional 
##  will first transform

p0.rel <- microbiome::transform(p0, "compositional")

grp_abund <- get_group_abundances(p0.rel, 
                                  level = "Phylum", 
                                  group="condition",
                                  transform = "compositional")

## An additonal column Sam_rep with sample names is created for reference purpose

mycols <- c("brown3", "steelblue","grey50")
# clean names 
grp_abund$OTUID <- gsub("p__", "",grp_abund$OTUID)
grp_abund$OTUID <- ifelse(grp_abund$OTUID == "", 
                          "Unclassified", grp_abund$OTUID)
mean.plot <- grp_abund %>% # input data
  ggplot(aes(x= reorder(OTUID, mean_abundance), # reroder based on mean abundance
             y= mean_abundance,
             fill=condition)) + # x and y axis 
  geom_bar(     stat = "identity", 
                position=position_dodge()) + 
  scale_fill_manual("condition", values=mycols) + # manually specify colors
  theme_bw() + # add a widely used ggplot2 theme
  ylab("Mean Relative Abundance") + # label y axis
  xlab("Phylum") + # label x axis
  coord_flip() # rotate plot 
mean.plot

## Plotting density of reads per sample

p0 <- microSubt
print_ps(p0)

## 01] ntaxa = 491

## 02] nsamples = 56

## 03] nsamplesvariables = 4

## 04] nranks = 6

## 05] Min. number of reads = 2766

## 06] Max. number of reads = 2766

## 07] Total number of reads = 154896

## 08] Average number of reads = 2766

## 09] Median number of reads = 2766

## 10] Sparsity = 0.792588012801862

## 11] Number of singletons = 179

## 12] % of taxa that are singletons 
##   (i.e. exactly one read detected across all samples) = 36.4562118126273

kable(head(tax_table(p0)))

	Kingdom	Phylum	Class	Order	Family	Genus
Otu086	Bacteria	Actinobacteriota	Actinobacteria	Bifidobacteriales	Bifidobacteriaceae	Bifidobacterium
Otu166	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	uncultured
Otu201	Bacteria	Actinobacteriota	Coriobacteriia	Coriobacteriales	Coriobacteriales_unclassified	Coriobacteriales
Otu170	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	GCA-900066575
Otu125	Bacteria	Firmicutes	Clostridia	Oscillospirales	Ruminococcaceae	Incertae
Otu512	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospiraceae

# reduce size for example
ps0 <- core(p0, detection = 10, prevalence = 20 / 100)
# Add a prefix to taxa labels
ps0.f2 <- format_to_besthit(ps0, prefix = "MyBug-")
kable(head(tax_table(ps0.f2))[3:6])

	Domain	Phylum	Class	Order	Family	Genus	best_hit
MyBug-Otu011:Lachnospiraceae	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospiraceae	MyBug-Otu011:Lachnospiraceae
MyBug-Otu010:Turicibacter	Bacteria	Firmicutes	Bacilli	Erysipelotrichales	Erysipelotrichaceae	Turicibacter	MyBug-Otu010:Turicibacter
MyBug-Otu038:Lachnospiraceae	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospiraceae	MyBug-Otu038:Lachnospiraceae
MyBug-Otu037:Oscillibacter	Bacteria	Firmicutes	Clostridia	Oscillospirales	Oscillospiraceae	Oscillibacter	MyBug-Otu037:Oscillibacter

p <- plot_read_distribution(ps0.f2, groups = "condition", 
                            plot.type = "density")
print(p + theme_biome_utils())

ps0 <- core(p0, detection = 10, prevalence = 20 / 100)
pseq_df <- phy_to_ldf(ps0, transform.counts = NULL)

## An additonal column Sam_rep with sample names is created for reference purpose

kable(head(pseq_df))

OTUID	Kingdom	Phylum	Class	Order	Family	Genus	Sam_rep	Abundance	group	site	mouse	condition
Otu046	Bacteria	Firmicutes	Clostridia	Oscillospirales	Ruminococcaceae	Ruminococcaceae	ce4406	34	ce4406	cecum	4406	10%Fat
Otu015	Bacteria	Firmicutes	Clostridia	Oscillospirales	Oscillospiraceae	Colidextribacter	ce4406	76	ce4406	cecum	4406	10%Fat
Otu011	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospiraceae	ce4406	132	ce4406	cecum	4406	10%Fat
Otu010	Bacteria	Firmicutes	Bacilli	Erysipelotrichales	Erysipelotrichaceae	Turicibacter	ce4406	0	ce4406	cecum	4406	10%Fat
Otu038	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Lachnospiraceae	ce4406	0	ce4406	cecum	4406	10%Fat
Otu037	Bacteria	Firmicutes	Clostridia	Oscillospirales	Oscillospiraceae	Oscillibacter	ce4406	27	ce4406	cecum	4406	10%Fat

# Plot Density of phyla
pseq <- microSubt

# check 10%Fat
control_ps <- subset_samples(pseq, condition=="10%Fat")
p_hc <- taxa_distribution(control_ps) + 
  theme_biome_utils() + 
  labs(title = "10%Fat")

# check Oleate
oleate_ps <- subset_samples(pseq, condition=="Oleate")
p_oleate <- taxa_distribution(oleate_ps) + 
  theme_biome_utils() + 
  labs(title = "Oleate")

# check 60%Fat
HFD_ps <- subset_samples(pseq, condition=="60%Fat")
p_hfd <- taxa_distribution(HFD_ps) + 
  theme_biome_utils() + 
  labs(title = "60%Fat")

# harnessing the power of patchwork
p_hc / p_oleate / p_hfd + plot_layout(guides = "collect") + 
  plot_annotation(tag_levels = "A")

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 210 rows containing non-finite values (stat_density).

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 184 rows containing non-finite values (stat_density).

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Removed 203 rows containing non-finite values (stat_density).

## Filtering the dataset

# Select specific taxa
microFirm <- subset_taxa(microSubt, Phylum %in% "Firmicutes")
random_tree2 = rtree(ntaxa(microFirm), rooted=TRUE, tip.label=taxa_names(microFirm))
plot(random_tree2)

microFirm = merge_phyloseq(microFirm, random_tree2)
microFirm

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 426 taxa and 56 samples ]
## sample_data() Sample Data:       [ 56 samples by 4 sample variables ]
## tax_table()   Taxonomy Table:    [ 426 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 426 tips and 1 internal nodes ]

# Subset samples by site or condition
microFecal <- subset_samples(microSubt, site=="fecal")
microFecal

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 491 taxa and 19 samples ]
## sample_data() Sample Data:       [ 19 samples by 4 sample variables ]
## tax_table()   Taxonomy Table:    [ 491 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 491 tips and 490 internal nodes ]

# if you needed to remove candidate outliers, can use the below to remove sample e.g. f4406
microSubt.1 <- prune_samples(sample_names(microSubt) != "f4406", microSubt)

Convert Phyloseq object to a dataframe

x<-psmelt(microSubt)
head(x)

##       OTU Sample Abundance group  site mouse condition  Kingdom     Phylum
## 11 Otu001  f4424      1571 f4424 fecal  4424    60%Fat Bacteria Firmicutes
## 50 Otu001  f4412      1220 f4412 fecal  4412    10%Fat Bacteria Firmicutes
## 24 Otu001  f4418      1196 f4418 fecal  4418    Oleate Bacteria Firmicutes
## 46 Otu001  f4407      1161 f4407 fecal  4407    10%Fat Bacteria Firmicutes
## 9  Otu001  f4406       933 f4406 fecal  4406    10%Fat Bacteria Firmicutes
## 53 Otu001  f4411       929 f4411 fecal  4411    10%Fat Bacteria Firmicutes
##      Class           Order           Family       Genus
## 11 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 50 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 24 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 46 Bacilli Lactobacillales Streptococcaceae Lactococcus
## 9  Bacilli Lactobacillales Streptococcaceae Lactococcus
## 53 Bacilli Lactobacillales Streptococcaceae Lactococcus

Session Info

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.1.1             ggpubr_0.4.0               
##  [3] RColorBrewer_1.1-2          viridis_0.5.1              
##  [5] viridisLite_0.3.0           microbiomeutilities_1.00.15
##  [7] dplyr_1.0.4                 microbiome_1.12.0          
##  [9] VennDiagram_1.6.20          futile.logger_1.4.3        
## [11] microbial_0.0.17            DESeq2_1.30.1              
## [13] SummarizedExperiment_1.20.0 Biobase_2.50.0             
## [15] MatrixGenerics_1.2.1        matrixStats_0.58.0         
## [17] GenomicRanges_1.42.0        GenomeInfoDb_1.26.2        
## [19] IRanges_2.24.1              S4Vectors_0.28.1           
## [21] BiocGenerics_0.36.0         ape_5.4-1                  
## [23] ellipse_0.4.2               clustsig_1.1               
## [25] knitr_1.31                  directlabels_2021.1.13     
## [27] vegan_2.5-7                 lattice_0.20-41            
## [29] permute_0.9-5               plyr_1.8.6                 
## [31] ggplot2_3.3.3               phyloseq_1.34.0            
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1             backports_1.2.1          fastmatch_1.1-0         
##   [4] igraph_1.2.6             splines_4.0.3            BiocParallel_1.24.1     
##   [7] digest_0.6.27            foreach_1.5.1            htmltools_0.5.1.1       
##  [10] fansi_0.4.2              magrittr_2.0.1           memoise_2.0.0           
##  [13] cluster_2.1.1            DECIPHER_2.18.1          openxlsx_4.2.3          
##  [16] limma_3.46.0             Biostrings_2.58.0        annotate_1.68.0         
##  [19] RcppParallel_5.0.3       prettyunits_1.1.1        jpeg_0.1-8.1            
##  [22] colorspace_2.0-0         ggrepel_0.9.1            blob_1.2.1              
##  [25] haven_2.3.1              xfun_0.21                crayon_1.4.1            
##  [28] RCurl_1.98-1.2           jsonlite_1.7.2           genefilter_1.72.1       
##  [31] dada2_1.18.0             phangorn_2.5.5           survival_3.2-7          
##  [34] iterators_1.0.13         glue_1.4.2               gtable_0.3.0            
##  [37] zlibbioc_1.36.0          XVector_0.30.0           DelayedArray_0.16.1     
##  [40] car_3.0-10               Rhdf5lib_1.12.1          abind_1.4-5             
##  [43] scales_1.1.1             pheatmap_1.0.12          futile.options_1.0.1    
##  [46] DBI_1.1.1                edgeR_3.32.1             rstatix_0.7.0           
##  [49] Rcpp_1.0.6               xtable_1.8-4             progress_1.2.2          
##  [52] foreign_0.8-81           bit_4.0.4                httr_1.4.2              
##  [55] ellipsis_0.3.1           farver_2.0.3             pkgconfig_2.0.3         
##  [58] XML_3.99-0.5             sass_0.3.1               locfit_1.5-9.4          
##  [61] utf8_1.1.4               labeling_0.4.2           tidyselect_1.1.0        
##  [64] rlang_0.4.10             reshape2_1.4.4           AnnotationDbi_1.52.0    
##  [67] munsell_0.5.0            cellranger_1.1.0         tools_4.0.3             
##  [70] cachem_1.0.4             cli_2.3.1                generics_0.1.0          
##  [73] RSQLite_2.2.3            ade4_1.7-16              broom_0.7.5             
##  [76] evaluate_0.14            biomformat_1.18.0        stringr_1.4.0           
##  [79] fastmap_1.1.0            yaml_2.2.1               bit64_4.0.5             
##  [82] zip_2.1.1                randomForest_4.6-14      purrr_0.3.4             
##  [85] nlme_3.1-152             formatR_1.7              rstudioapi_0.13         
##  [88] compiler_4.0.3           curl_4.3                 png_0.1-7               
##  [91] ggsignif_0.6.1           tibble_3.0.6             geneplotter_1.68.0      
##  [94] bslib_0.2.4              stringi_1.5.3            highr_0.8               
##  [97] forcats_0.5.1            Matrix_1.3-2             multtest_2.46.0         
## [100] vctrs_0.3.6              pillar_1.5.0             lifecycle_1.0.0         
## [103] rhdf5filters_1.2.0       jquerylib_0.1.3          data.table_1.14.0       
## [106] bitops_1.0-6             R6_2.5.0                 latticeExtra_0.6-29     
## [109] hwriter_1.3.2            ShortRead_1.48.0         gridExtra_2.3           
## [112] rio_0.5.16               codetools_0.2-18         lambda.r_1.2.4          
## [115] MASS_7.3-53.1            assertthat_0.2.1         rhdf5_2.34.0            
## [118] withr_2.4.1              GenomicAlignments_1.26.0 Rsamtools_2.6.0         
## [121] GenomeInfoDbData_1.2.4   mgcv_1.8-34              hms_1.0.0               
## [124] gghalves_0.1.1           quadprog_1.5-8           tidyr_1.1.2             
## [127] rmarkdown_2.7            carData_3.0-4            Rtsne_0.15